function breakpoints = find_breakpoints(img)

% threshold image and find edges
img = ~im2bw(img, 0.5);

% find first non-zero pixel
[y x] = find(img);
starty = y(1);
startx = x(1);
P = [starty     startx];

% trace boundary
B = bwtraceboundary(img, P, 'E');

% setup parameters
k2 = 3;
k3 = 4;
k4 = 0.6;

% loop over contour and find curvature
curvature = zeros(size(B, 1), 1);

for i=1:size(B, 1)
    a = [B(mod(i + k2 - 1, size(B, 1)) + 1, 2) - B(i, 2); B(mod(i + k2 - 1, size(B, 1)) + 1, 1) - B(i, 1)];
    b = [B(mod(i - k2 - 1, size(B, 1)) + 1, 2) - B(i, 2); B(mod(i - k2 - 1, size(B, 1)) + 1, 1) - B(i, 1)];
    curvature(i) = (a' * b) / (norm(a) * norm(b));
end

% find apparent breakpoints

breakpoints = [];
for i=1:size(curvature, 1)
    first = mod(i - k3 -1, size(curvature, 1)) + 1;
    last = mod(i + k3 - 1, size(curvature, 1)) + 1;
    fun = wrapIndex(curvature, first, last);
    
    if (curvature(i) == max(fun) && curvature(i) > -k4)
        breakpoints = [breakpoints i];
    end
end

% convert breakpoint indices to breakpoint locations
breakpoints = B(breakpoints', :);

% show breakpoints and image
imshow(img);
hold on;
plot(breakpoints(:, 2)', breakpoints(:, 1)', 'b+', 'MarkerSize', 15);

end

function res = wrapIndex(mat, first, last)
    if first < last
        res = mat(first:last);
    else
        res = [mat(first:end) ; mat(1:last)];
    end
end